
do 99_programs.do

** continuous fines, QQ desired at different fine percentiles
** generate dta for Figure A4 and A5
use ../data/c8290hh_reg_main, clear

* fine polynomials
g fine3b2 = fine3b * fine3b
g fine3b3 = fine3b2 * fine3b
g fine3b4 = fine3b3 * fine3b
g fine3b5 = fine3b4 * fine3b

g fine3b2_twinhh = fine3b2 * twinhh
g fine3b3_twinhh = fine3b3 * twinhh
g fine3b4_twinhh = fine3b4 * twinhh
g fine3b5_twinhh = fine3b5 * twinhh


local ctrl1 = " age2birthP2 male age age2 i.edulevP2 i.prov i.birthyP2 i.wave prov*_birthyP2 "

* poly 3
est clear

local ctrl = " `ctrl1'  "
eststo: bootstrap _b beta_diff_20 = (_b[beta_B_20] - _b[beta_A]) beta_diff_30 = (_b[beta_B_30] - _b[beta_A]) ///
beta_diff_40 = (_b[beta_B_40] - _b[beta_A]) beta_diff_50 = (_b[beta_B_50] - _b[beta_A]) ///
beta_diff_60 = (_b[beta_B_60] - _b[beta_A]) beta_diff_70 = (_b[beta_B_70] - _b[beta_A]) ///
beta_diff_80 = (_b[beta_B_80] - _b[beta_A]) beta_diff_90 = (_b[beta_B_90] - _b[beta_A]) ///
beta_diff_100 = (_b[beta_B_100] - _b[beta_A]) beta_diff_110 = (_b[beta_B_110] - _b[beta_A]) ///
beta_diff_120 = (_b[beta_B_120] - _b[beta_A]), ///
seed(123) rep(100) cluster(prov_edulevP2) strata(prov) : reg_nchild3_edulev3_fine3b3 `ctrl'


local kvars1 "alpha_z alpha_zx1 alpha_zx2 alpha_zx3 alpha_x1 alpha_x2 alpha_x3 " 
local kvars2 "rho_z rho_zx1 rho_zx2 rho_zx3 rho_x1 rho_x2 rho_x3  beta_A r2_fs r2_rf "
local kvars3 "beta_B_20 beta_B_30 beta_B_40 beta_B_50 beta_B_60 beta_B_70 "
local kvars4 "beta_B_80 beta_B_90 beta_B_100 beta_B_110 beta_B_120 "
local kvars5 "beta_diff_20 beta_diff_30 beta_diff_40 beta_diff_50 beta_diff_60 beta_diff_70 "
local kvars6 "beta_diff_80 beta_diff_90 beta_diff_100 beta_diff_110 beta_diff_120 "

local kvars "`kvars1' `kvars2' `kvars3' `kvars4' `kvars5' `kvars6' "
regsave `kvars' using ../Rdta/finepoly3, ci level(90) replace


* poly 4
est clear

local ctrl = " `ctrl1'  "
eststo: bootstrap _b beta_diff_20 = (_b[beta_B_20] - _b[beta_A]) beta_diff_30 = (_b[beta_B_30] - _b[beta_A]) ///
beta_diff_40 = (_b[beta_B_40] - _b[beta_A]) beta_diff_50 = (_b[beta_B_50] - _b[beta_A]) ///
beta_diff_60 = (_b[beta_B_60] - _b[beta_A]) beta_diff_70 = (_b[beta_B_70] - _b[beta_A]) ///
beta_diff_80 = (_b[beta_B_80] - _b[beta_A]) beta_diff_90 = (_b[beta_B_90] - _b[beta_A]) ///
beta_diff_100 = (_b[beta_B_100] - _b[beta_A]) beta_diff_110 = (_b[beta_B_110] - _b[beta_A]) ///
beta_diff_120 = (_b[beta_B_120] - _b[beta_A]), ///
seed(123) rep(100) cluster(prov_edulevP2) strata(prov) : reg_nchild3_edulev3_fine3b4 `ctrl'


local kvars1 "alpha_z alpha_zx1 alpha_zx2 alpha_zx3 alpha_zx4 alpha_x1 alpha_x2 alpha_x3 alpha_x4 " 
local kvars2 "rho_z rho_zx1 rho_zx2 rho_zx3 rho_zx4 rho_x1 rho_x2 rho_x3 rho_x4 beta_A r2_fs r2_rf "
local kvars3 "beta_B_20 beta_B_30 beta_B_40 beta_B_50 beta_B_60 beta_B_70 "
local kvars4 "beta_B_80 beta_B_90 beta_B_100 beta_B_110 beta_B_120 "
local kvars5 "beta_diff_20 beta_diff_30 beta_diff_40 beta_diff_50 beta_diff_60 beta_diff_70 "
local kvars6 "beta_diff_80 beta_diff_90 beta_diff_100 beta_diff_110 beta_diff_120 "

local kvars "`kvars1' `kvars2' `kvars3' `kvars4' `kvars5' `kvars6' "
regsave `kvars' using ../Rdta/finepoly4, ci level(90) replace


* poly 5
est clear

local ctrl = " `ctrl1'  "
eststo: bootstrap _b beta_diff_20 = (_b[beta_B_20] - _b[beta_A]) beta_diff_30 = (_b[beta_B_30] - _b[beta_A]) ///
beta_diff_40 = (_b[beta_B_40] - _b[beta_A]) beta_diff_50 = (_b[beta_B_50] - _b[beta_A]) ///
beta_diff_60 = (_b[beta_B_60] - _b[beta_A]) beta_diff_70 = (_b[beta_B_70] - _b[beta_A]) ///
beta_diff_80 = (_b[beta_B_80] - _b[beta_A]) beta_diff_90 = (_b[beta_B_90] - _b[beta_A]) ///
beta_diff_100 = (_b[beta_B_100] - _b[beta_A]) beta_diff_110 = (_b[beta_B_110] - _b[beta_A]) ///
beta_diff_120 = (_b[beta_B_120] - _b[beta_A]), ///
seed(123) rep(100) cluster(prov_edulevP2) strata(prov) : reg_nchild3_edulev3_fine3b5 `ctrl'


local kvars1 "alpha_z alpha_zx1 alpha_zx2 alpha_zx3 alpha_zx4 alpha_zx5 alpha_x1 alpha_x2 alpha_x3 alpha_x4 alpha_x5 " 
local kvars2 "rho_z rho_zx1 rho_zx2 rho_zx3 rho_zx4 rho_zx5 rho_x1 rho_x2 rho_x3 rho_x4 rho_x5  beta_A r2_fs r2_rf "
local kvars3 "beta_B_20 beta_B_30 beta_B_40 beta_B_50 beta_B_60 beta_B_70 "
local kvars4 "beta_B_80 beta_B_90 beta_B_100 beta_B_110 beta_B_120 "
local kvars5 "beta_diff_20 beta_diff_30 beta_diff_40 beta_diff_50 beta_diff_60 beta_diff_70 "
local kvars6 "beta_diff_80 beta_diff_90 beta_diff_100 beta_diff_110 beta_diff_120 "

local kvars "`kvars1' `kvars2' `kvars3' `kvars4' `kvars5' `kvars6' "
regsave `kvars' using ../Rdta/finepoly5, ci level(90) replace




*** Table A7: gender-interactive model
use ../data/c8290hh_reg_main, clear

g byte female = 1 - male

est clear

local ctrl1_nomale = " age2birthP2   age age2 i.edulevP2 i.prov i.birthyP2 i.wave prov*_birthyP2 "

local ctrl = " `ctrl1_nomale' "
eststo: bootstrap _b beta_diff = (_b[beta_B] - _b[beta_A]) beta_diff_f = (_b[beta_B_f] - _b[beta_A_f]) ///
beta_A_gender_dif = (_b[beta_A] - _b[beta_A_f]) beta_B_gender_dif = (_b[beta_B] - _b[beta_B_f]) ///
beta_did = ( (_b[beta_B] - _b[beta_A]) - (_b[beta_B_f] - _b[beta_A_f]) ), ///
seed(123) rep(100) cluster(prov_edulevP2) strata(prov) : reg_nchild3_edulev3_gender `ctrl'


local kvar1 "alpha_z alpha_zx alpha_x alpha_z_f alpha_zx_f alpha_x_f alpha_f "
local kvar2 "rho_z rho_zx rho_x rho_z_f rho_zx_f rho_x_f rho_f "
local kvar3 "beta_A beta_B beta_A_f beta_B_f beta_diff beta_diff_f beta_A_gender_dif beta_B_gender_dif beta_did r2_fs r2_rf "

local kvars "`kvar1' `kvar2' `kvar3' "
esttab using ../output/TableA7_tabGenderInteract.csv, ///
keep(`kvars')  /// 
order(`kvars') ///
b(3) se(3) star(* .10 ** .05 *** .01) stats(N, labels("Observations") fmt(0)) ///
title("Rural second-born twins") nogaps label n replace



***** Table A8: higher parities


*** nchild4, sample3plus
use ../data/c8290hh, clear

keep if birthyP2 >= 1940 & birthyP2 <= 1960
keep if age >= 13 & age <= 17
keep if sample3plus == 1
keep if urban == 0
keep if han == 1

drop if age1birthP2 <= 15

sum fine4b, d
drop if fine4b > r(p99)

g fine4b_twinhh = fine4b * twinhh

g male_twinhh = male * twinhh

g eduyP2dm = eduyP2 - 6
g eduyP2dm_twinhh = eduyP2dm * twinhh

g twinsex = 1 if nmaletwin == 2 & twinhh
replace twinsex = 2 if nmaletwin == 0 & twinhh
replace twinsex = 3 if nmaletwin == 1 & twinhh

tab twinsex

* province-specific trends, nnumber of distinct prov - 1
tab prov, g(prov)
distinct prov
local n = `r(ndistinct)' - 1
forval i=1/`n' {
	g prov`i'_birthyP2 = prov`i' * birthyP2
	g prov`i'_year2 = prov`i' * (birthyP2 - 1940) * (birthyP2 - 1940)
}

egen prov_edulevP2 = group(prov edulevP2)

	* log
foreach var of varlist lossratio CRdratio sdyReceivedRatio minShare gdppc popGrowth {
	gen ln`var' = log(`var')
}

	* de-median fine3b and baseline controls
foreach var of varlist fine4b age2birthP2 age age2 eduyP2 {
	sum `var', d
	g `var'_dm = `var' - r(p50)
}

tab twinorder

	* summary statistics
tabstat edulev3 nchild4 fine4b twinhh  male age  eduyP2 ageP2 age1birthP2 age3birthP2, s(mean sd n) by(twinhh)


local ctrl1_4 = " age3birthP2 male age age2 i.edulevP2 i.prov i.birthyP2 i.wave prov*_birthyP2 "

est clear

local ctrl = " `ctrl1_4'  "
eststo: bootstrap _b beta_diff = (_b[beta_B] - _b[beta_A]), ///
seed(123) rep(100) cluster(prov_edulevP2) strata(prov) : reg_nchild4_edulev3 `ctrl'

local kvars "alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B alpha_sum beta_diff r2_fs r2_rf "
esttab using ../output/TableA8_col1.csv, ///
keep(`kvars')  /// 
order(`kvars') ///
b(3) se(3) star(* .10 ** .05 *** .01) stats(N, labels("Observations") fmt(0)) ///
title("Rural second-born twins") nogaps label n replace

est clear

local ctrl = " `ctrl1_4'  "
eststo: bootstrap, seed(123) rep(100) cluster(prov_edulevP2) strata(prov): reg_nchild4 `ctrl'
test twinhh fine4b_twinhh, df(104)
estadd scalar F = r(F)

local kvars "twinhh fine4b_twinhh fine4b "
esttab using ../output/TableA8_col1_F.csv, ///
keep(`kvars')  /// 
order(`kvars') ///
b(3) se(3) star(* .10 ** .05 *** .01) stats(F N, labels("F-statistic" "Observations") fmt(2 0)) ///
title("Rural second-born twins") nogaps label n replace





*** nchild5, sample4plus
use ../data/c8290hh, clear

keep if birthyP2 >= 1940 & birthyP2 <= 1960
keep if age >= 13 & age <= 17
keep if sample4plus == 1
keep if urban == 0
keep if han == 1

drop if age1birthP2 <= 15

sum fine5b, d
drop if fine5b > r(p99)

g fine5b_twinhh = fine5b * twinhh

g male_twinhh = male * twinhh

g eduyP2dm = eduyP2 - 6
g eduyP2dm_twinhh = eduyP2dm * twinhh

g twinsex = 1 if nmaletwin == 2 & twinhh
replace twinsex = 2 if nmaletwin == 0 & twinhh
replace twinsex = 3 if nmaletwin == 1 & twinhh

tab twinsex


* province-specific trends, nnumber of distinct prov - 1
tab prov, g(prov)
distinct prov
local n = `r(ndistinct)' - 1
forval i=1/`n' {
	g prov`i'_birthyP2 = prov`i' * birthyP2
	g prov`i'_year2 = prov`i' * (birthyP2 - 1940) * (birthyP2 - 1940)
}

egen prov_edulevP2 = group(prov edulevP2)

	* log
foreach var of varlist lossratio CRdratio sdyReceivedRatio minShare gdppc popGrowth {
	gen ln`var' = log(`var')
}

	* de-median fine3b and baseline controls
foreach var of varlist fine5b age2birthP2 age age2 eduyP2 {
	sum `var', d
	g `var'_dm = `var' - r(p50)
}

tab twinorder

local ctrl1_5 = " age4birthP2 male age age2 i.edulevP2 i.prov i.birthyP2 i.wave prov*_birthyP2 "

est clear

local ctrl = " `ctrl1_5'  "
eststo: bootstrap _b beta_diff = (_b[beta_B] - _b[beta_A]), ///
seed(123) rep(100) cluster(prov_edulevP2) strata(prov) : reg_nchild5_edulev3 `ctrl'

local kvars "alpha_z alpha_zx alpha_x rho_z rho_zx rho_x beta_A beta_B alpha_sum beta_diff r2_fs r2_rf "
esttab using ../output/TableA8_col2.csv, ///
keep(`kvars')  /// 
order(`kvars') ///
b(3) se(3) star(* .10 ** .05 *** .01) stats(N, labels("Observations") fmt(0)) ///
title("Rural second-born twins") nogaps label n replace

est clear

local ctrl = " `ctrl1_5'  "
eststo: bootstrap, seed(123) rep(100) cluster(prov_edulevP2) strata(prov): reg_nchild5 `ctrl'
test twinhh fine5b_twinhh, df(104)
estadd scalar F = r(F)

local kvars "twinhh fine5b_twinhh fine5b "
esttab using ../output/TableA8_col2_F.csv, ///
keep(`kvars')  /// 
order(`kvars') ///
b(3) se(3) star(* .10 ** .05 *** .01) stats(F N, labels("F-statistic" "Observations") fmt(2 0)) ///
title("Rural second-born twins") nogaps label n replace



